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Abstract 

DNA nanocompartmcnt is a typical DNA-based machine whose func- 
tion is dependent of molecular collective effect. Fundamental prop- 
erties of the device have been addressed via electrochemical analysis, 
fluorescent microscopy, and atomic force microscopy. Interesting and 
novel phenomena emerged during the switching of the device. We have 
found that DNAs in this system exhibit a much steep melting transi- 
tion compared to ones in bulk solution or conventional DNA array. 
To achieve an understanding to this discrepancy, we introduced DNA- 
DNA interaction potential to the conventional Ising-like Zimm-Bragg 
theory and Peyrard-Bishop model of DNA melting. To avoid unrealis- 
tic numerical calculation caused by modification of the Peyrard-Bishop 
nonlinear Hamiltonian with the DNA-DNA interaction, we established 
coarse-gained Monte Carlo recursion relations by elucidation of five 
components of energy change during melting transition. The result 
suggests that DNA-DNA interaction potential accounts for the ob- 
served steep transition. 



1 .Introduction 

Studies on the physical chemistry of DNA denaturation have been lasted 
for almost forty years [1-3]. In 1964, Lifson proposed that a phase tran- 
sition exists in one-dimensional polymer structure. He introduced several 
pivotal concepts, like sequence partition function, sequence generating func- 
tion, etc., and established a systematic method to calculate the partition 
function [1]. These allow us to derive important thermodynamic quantities 
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of the system. In 1966, Poland and Scherage applied Lifson's method to 
conduct research on amino acid and nucleic acid chains. They built Poland- 
Scherage (PS) model for calculating the sequence partition function and 
discussing the behavior of polymers in melting transitions. 

Another excellent progress would be the building of Peyrard-Bishop (PB) 
model [4,5] for DNA chains. In PB model, the Hamiltonian of a single DNA 
chain, which is constructed by phonon calculations, is given so that we can 
obtain the system properties through statistical physics method. The PB 
model has introduced mathematical formula of stacking energy, as well as 
the kinetic energy and potential energy of each base pair. By theoretical 
calculation, one can show the entropy-driven transition that leads DNA to 
shift from ordered state to disorder one [6,7]. 

However, all these works have not involved the DNA-DNA interactions 
because the subject investigated is DNAs in bulk solution, and the interac- 
tion between them has ever been neglected. The main idea of this paper 
is to inspect the influence of collective effect on the DNA melting process, 
primarily motivated by the experiment results of DNA nanocompartment 
[8,9]. Under the enlightenment of Poland- Scherage model and Zimm-Bragg 
model [10], we simplify Peyrard-Bishop model to meet a reasonable Monte 
Carlo simulation by the elucidation of five components of energy changes 
during melting transition. The result shows that the melting temperature 
and transition duration depend on whether we take into account the DNA- 
DNA interactions among columnar assemblies of DNA. 

2. Experiment 

Recently, we found that specially designed DNA array can form a molecular 
cage on surfaces [8,9]. This molecular cage is switchable due to allosteric 
transformation driven by the collective hybridization of DNA. We named 
it "active DNA nanocompartment (ADNC)". Typical DNA motif designed 
to fabricate ADNC comprises two contiguous elements (inset to figure la): 
a double-stranded DNA (dsDNA) whose array is responsible for a compact 
membrane (figure la, right), and a single-stranded DNA (ssDNA) serving as 
skeleton supporting the dsDNA membrane, which is terminated on its 5 end 
by a surface linker such as an alkanethiol group that can be tethered to gold 
surface with a sulphur-gold bond [9] or an amino group that can be tethered 
to Si02 substrate with specific surface attachment chemistry [11]. Because 
the diameter of ssDNA is much smaller than that of dsDNA, a compart- 
ment with designable effective height (heff, 5 ~ 50nm, commensurate with 
the length of ssDNA skeleton) can form between the dsDNA membrane and 
substrate surface. 

Since ADNC is reversibly switchable, it is able to encage molecules with 
suitable size. We name this phenomenon molecular encaging effect. Both 
electrochemical methods [12] and fluorescent microscopy are used to sub- 
stantiate the molecular encaging effect and the reversibility of switching. 
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Once the closed ADNC entraps some chemical reporters, the surface concen- 
tration (r nc ) of the encaged reporters can be determined by cyclic voltammc- 
try or fluorescent microscopy. Figure lb shows the isotherms of the molec- 
ular encaging effect for fluorescein (C2oHioNa20s). Figure lc presents the 
melting curves of ADNC. Using the encaged molecules as indicator greatly 
sharpens the melting profiles for the perfectly complementary targets, and 
flattens denaturation profiles for the strands with a wobble mismatch. The 
observation shows that single-base mismatched strands are incapable of clos- 
ing ADNC on surfaces. The result is highly consistent to our observation 
by electrochemical analysis [12]. These observations bring up an intriguing 
question: why the melting curves exhibit so steep transition compared to 
the case of DNA in bulk solutions or on a loosely packed microarray? We 
try to address this question in this paper. 

Worthy of mention is that the steepness of melting transition is use- 
ful when the ADNC is applied to DNA detection [8,9]. First, it greatly 
enhances the discrepancy of perfect targets and single mismatches. This 
provides much enhanced specificity in DNA recognition, 100 : 1 ~ 105 : 1 of 
our system versus 2.7 : 1 of conventional system. Second, more sensitivity is 
obtained with optimally decreased ambiguity. Therefore, the clarification of 
the origin of the steep shape should help us to further extend the experience 
to related fields or generate new techniques. 

3.Modeling 

Taking into account the directional specificity of the hydrogen bonds, the 
Hamiltonian of a single DNA chain is obtained as following form according 
to PB model [4-6], 



^my n 2 + w{y n ,y n ^i) + V(y„ 



(1) 



where the y n is the component of the relative displacement of bases along the 
direction of hydrogen bond. The stacking energy w(y n ,y n -i) corresponds 
to the interaction between neighboring base pair in one DNA chain 



k 

w(y n ,y n -i) = g 



1 + pe 



-a(y„+y„-i) 



(Vn ~ Vn-l? (2) 



The Morse potential describes the potential for the hydrogen bonds 

V = D{e~ ay - l) 2 (3) 

However, in this study, the Hamiltonian in equation (1) is not sufficient; 
it neglects the structure of close-packing of DNA in ADNC. In our system, 
one should take into account the interactions between the nearest neighbor- 
ing molecules [13,14]. To model the interaction, one envisions the molecules 
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Figure 1: (a) Schematic drawing of a dual-state ADNC. By adding or re- 
moving 'fuel' strands (shorter segment), the ADNC can be switched between 
on (right) and off (left) state. Inset, a typical sequence used to fabricate 
ADNC. The 'fuel' strands is a segment of human p53 gene containing one 
site of most frequent mutation, (b) Isotherms of T for fluorescein encaged 
in a closed ADNC (h e ff = 15bp). The isotherm fits well to the Langmuir 
model: x/T nc = (1/T nc , max )x + (l/KT nCimax ), where x is the concentra- 
tion of the reporter and K is the association constant per site, T nc the 
surface concentration of encaged molecules, (c) Melting curves using the 
encaged fluorescein molecules as indicators. Filled circles or squares are cor- 
responding to perfect complementary strands, and hollow circles or squares 
to single-base mismatched strands. The unit of relative intensity of fluo- 
rescein is defined as the light intensity of 5/iZ WnM fluorescein on a spot 
size with 5mm diameter. Inset shows the length of the ssDNA skeleton of 
nanocompartment and the complementary type. The background noise is 
within 1 unit. 
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as rigid cylinders, carrying helical and continuous line charges on their sur- 
faces. Each DNA duplex carries the negative charge of phosphates plus a 
compensating positive charge from the adsorbed counterions. Let < 9 < 1 
be the degree of charge compensation, /i, h and h the fractions of con- 
densed counterions in the minor and major grooves (fi + h + h = 1). The 
mobile counterions in solution screen the Coulomb interactions between the 
two molecules, causing at large separations an exponential decay of the lat- 
ter with the Debye screening length re -1 . The solvent is accounted for by 
its dielectric constant e. The structural parameters of B-DNA are half az- 
imuthal width of the minor groove 4> s « 0.4vr , pitch H ss Mk(g = 2ir/H), 
and hard-core radius a = 9A. We take the following form for the pair 
interaction potential [15-18]: 

he + (-i) n f 2 e - (i - he) cos(n&)l 2 



u(R,4>) = u 

n=—c 

(-1)" cos(ngAz)K (K n R) - Q n ^ n (n n R, n n a 
X {K n /Kf[K' n {K n a)Y 



(4) 



where R(> 2a) is the distance between the two parallel DNA molecules, 
Az a vertical displacement, equivalent to a "spin angle" (p = gz. Here, 
U q = 8ira 2 /eK, 2 (about 2.9fc B T/A at physiological ionic strength), and n n — 
^k 2 + n 2 g 2 . n„, m (x, y) is given by 



^"n,m{x,y) — ^ ' 



I'M 



(5) 



with the modified Bessel functions K n {x) and Ijiy)- The primes denote 
derivatives. The sum rapidly converges, and it can be truncated after 
|n| = 2. Since K n R > 3 and g ~ k, each of the terms in the sum de- 
creases exponentially at increasing R with the decay length k^ 1 oc 1/n. 

Figure 2 present a scheme of interaction between two neighboring colum- 
nar DNA molecules charged with counterions on its surface. The distance 
between two DNA columns in our simulation is about 30A and the helical 
pitch of DNA molecule is about 36 A. For brevity, we take the mean-field 
approximation that the pair interactions mainly exist between charges in 
the same height. 



4. Monte Carlo Simulation 

Let t be the dimensionless variable to mark the time series of simulation 
(t = 0, 1, 2 . . .) and T the environmental temperature. Assuming that M xJV 
DNAs are on the ADNC, the position of each DNA can be represented by 
its coordinates (x, y), where x,y G N, and 0<x<M, < y < N. All DNA 
molecules in ADNC have identical sequence with P base pairs. Therefore 
there is the collection of MxiVxP base pairs. The degree of freedom of the 
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Figure 2: The pair interaction between two parallel B-DNA double helixes. 
The black balls with positive signs in the center represent the absorbed 
positive counterions, while the little grey balls represent the phosphate car- 
rying negative charges. Each DNA duplex carries the negative charge of 
phosphates with area density of 16.8fiC/cm 2 plus a compensating positive 
charge coming from the adsorbed counterions. We take the assumption that, 
and the distance between them are approximately 30A. 
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system is also M x N x P. The position of each base pair is thus represented 
by coordinates (x, y, i), where % € N, < i < P. We take that the indices of 
base pairs is assigned from the bottom to the top of the DNA. 

At the time to, the state for an arbitrary base pair at (xo,yo,i) with 
well-formed hydrogen bonds is represented as (xo,yo,i,to) = 1. Contrar- 
ily, the state of a base pair with decoupled hydrogen bonds is denoted as 
(xo,yo,i,to) = [10, 19]. ip(x,y,i,t) is a function of the time and the posi- 
tion of the base pair. Therefore, the state of each DNA molecules in ADNC 
can be represented by a sequence of digits. The number of all possible states 
is 2 MxNxP . 




11 1010 1000 1101 1100 0100 110 



Figure 3: A schematic graph for a double-stranded DNA molecule associ- 
ated with bool quantities to represent well-formed hydrogen bonds (denoted 
as 1) and decoupled bonds (denoted as 0). 



The simulation begins ati = 0,T = 0°C. At each step, t increases 
by 1, and state of base pair at (xo,yo,i) is inverted, i.e. ip(xo,yo,i,t + 
1) = 1 — (xo,yo,i,t). We assume that by changing the state of the system 
for M x N x P x Z times, the system will approximate the equilibrium 
state infinitely. The change will be applied to each base pair for average 
Z times. Z is determined by experience and should be reasonable. We 
increase T by AT during the simulation. Therefore we have the relation 
T = t AT 

Whether the state inversion is permitted depends on the energy change 
(AE t ) in each step. The possibility of the state change at each step is 



P^(t)^^(t + !)) = { 



1 for AE t < 

(6) 

e -AE t /k B T j Qr AEf > 



If the current state change is permitted, we keep up changing the system 
state at t + 1. If the state change is forbidden by the possibility, the system 
state remains unchanged at t and waits for another change at t + 1. 

To achieve a relatively precise simulation, the change of the total energy 
at time t + 1 relative to that at the time t is analyzed by five components. 
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The recursion relation of energy change in each step is written as: 



5 

AE t = E(t + 1) - E(t) = A£j (7) 

i=i 

where E(t+1) and E(t) are the system energy for the instant t + 1 and t 
respectively, and AEi (I = 1, 2, 3, 4, 5) is the variation of the Ith component. 

The energy change depends on both the recursion relation of the base 
pair at (xo,yo,i) and the states of its nearest neighbors. The global energy 
variation is determined by the local states around the base pair (xo,yo,i). 
Following analysis presents the recursion relation of energy changes. 

4-1. The hydrogen-binding energy (AE\) 

This component of the energy consists of the Morse potential (equation (3)) 
and kinetic energy along the orientation of the hydrogen bonds. The binding 
energy is independent of the states of its neighboring base pairs. 



AE l = < 



J if i>(x ,y ,i,t) = 

(8) 

-J if t/j(x ,y ,i,t) = 1 



where J(J < 0) is the binding energy for each base pair that is in '1' state, 
while the binding energy for the '0' state is zero to be reference. 

4-2. The stacking energy (AE2) 

To simplify the calculation of the stacking energy shown in equation (2), we 
take into account the states of base pairs at (xo, yo,i — l,t) and (xo, yo,i + 
1, t). Their states remain unchanged during the interval from ttot + 1. We 
employ the periodic boundary condition (PBC) listed below. 

ip(x + M, y, i, t) = ip(x, y, i, t) 

ip(x,y + N,i,t) = ip(x,y,i,t) (9) 
ip(x,y,i + P,t) = ip(x,y,i,t) 

Therefore i[)(xo,yo, i — l,t) and i/j( x o, Uo, i + l,t) are both well defined. The 
stacking energy reflects the interaction between nearest neighboring base 
pairs in same DNA, and it exists only when two nearest neighbors are in 
'1' state at the same time. We use the symbol ip(x, y, i\, i2, ■ ■ ■ , i n , t) = 
bi, 62, . . . , b n to denote states in the same DNA for convenience, which means 
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(x,y,h,t) = 61, (x, 
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(10) 



where w is the stacking energy stored in two nearest neighboring base pairs 
in '1' state. 



4-3. Morse potential away from equilibrium point (AE3) 
We set AEi = for the uncoupled hydrogen bond at base pairs. However, 
for a '0' state is next near to a '1' state in the same DNA strand, the distance 
between two base pairs is so close that the Morse potential should be taken 
into account. We assigned energy E to every two nearest neighboring base 
pairs that are in different states in the same DNA. 
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4-4- The effect of excluded volume (AE4) 

The effect of excluded volume in the nature of DNA phase transition is dis- 
cussed in Fisher's work [20]. The excluded volume effect is connected to the 
system entropy variation. The effect is prone to separate two complemen- 
tary strands in a double helix. We use F to represent the energy change 
corresponding to this effect. One should notice ^ < 0. We then have 
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The energy changes discussed above are summarized in the figure 4 be- 
low, which does not take into account the DNA-DNA interactions so far. 




Figure 4: Scheme for the energy changes. We calculate the energy change 
under every case when the state of a base pair takes inversion. The value 
overhanging the double-arrow represents the energy variation when the state 
changes from left to right A = AE X + AE 2 + AE 3 + AE 4 . 



4.5. DNA-DNA interaction potential (AE 5 ) 

We have introduced DNA-DNA interaction in previous section. For each 
base pair, we denote the state of its m nearest neighbors with Aj, (i = 
1, 2, . . . , m; Aj = 0, 1). AE§ can be written as 



AE, 



GTT=ih if *P(x ,y ,i,t) = 

(13) 

-GE™i if Y>(x ,yo,M) = 1 



where G is the interaction energy between each pair of ions. Adding AE5 
to A, we will get the energy variation including the DNA-DNA interaction. 

5. Results and Discussions 

Following the equation (5) - (12), we could achieve a coarse-gained simula- 
tion of the melting curves of ADNC as well as that of DNA in bulk solutions. 
To perform the task, we choose suitable scale parameters to carry out the 
simulation: M = 100, N = 100, P = 20. The values of M and N chosen 
are much smaller than ones of the actual situation, which is up to 10 4 in the 
experiment. Since we take the periodic boundary condition, the values of M 
and N used do not change our result. The starting temperature is 0°C, and 
the final temperature is 100 °C, with increment of 0.01 °C for each step. To 
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guarantee the system reaches equilibrium state, we take state changes under 
a specific temperature. Each base pair has average 5 times to be changed. 
At each step, we count the number of DNAs that is still hybridized and 
calculate the percentage for dsDNA in ADNC. The simulation result shown 
in figure 5 shows a steep melting transition (hollow circles), consistent to 
the experimental observations. The simulated result without considering 
DNA-DNA interaction show in filled circles in figure 5 also agrees with the 
DNA melting curves in bulk solution. Comparison between the two cases 
suggests that the DNA-DNA interaction greatly increases the melting point 
of dsDNA chains. 




20 40 60 80 100 

TemperatureCC) 



Figure 5: Simulation results of melting curves for collective DNA molecules. 
The filled circle represents the phase transition curve without considering 
the DNA-DNA interaction, while the hollow circle is the counterpart that 
takes into account the interaction. Parameters used in the simulations are: 
J = -1900&B, w = -250k B , E = 850k B , F = -1650^, = -10k B , 
G = -125k B . 

In conclusion, we have established a simple coarse-gained model to sim- 
ulate the melting transition of DNA in ADNC. The result provides a reason- 
able explanation for our experimental observations. Although the simula- 
tion method discretizes the Morse potential and stacking energy proposed in 
Peyrard-Bishop model, the result still present a comparable approximation 
to experimental data due to our fine treatment of energy changes during 
melting transition. However, this work is only the beginning of insight- 
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ful theoretical investigation for the rationality of ADNC. In future work, 
we will establish a more precise model to employ an extensive investigation 
of phase transition occurring in ADNC as well as its derived DNA machines. 
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